Fig 4 Postmortem

Volcano

nygc_postmortem_spinal_cord.als.labels <- nygc_postmortem_spinal_cord$res %>% filter(gene_name %in% ALS_genes) %>% top_n(n=5,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.p53.labels <- nygc_postmortem_spinal_cord$res %>% filter(gene_name %in% p53.DDR.repair_gene_set) %>% top_n(n=5,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.labels <- nygc_postmortem_spinal_cord$res %>% filter(!grepl("^AC|^AL", gene_name)) %>% top_n(n=10,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.volcano <- volcano_plot(nygc_postmortem_spinal_cord$res, numerator = "ALS", denomenator = "CTRL",
               title = "Postmortem spinal cord", subtitle = "57 Controls vs 214 ALS", ymax = 45, xmax = 2.9, dpi = 300,
               gene_labels = c(nygc_postmortem_spinal_cord.als.labels, #nygc_postmortem_spinal_cord.ipsn.up.labels, nygc_postmortem_spinal_cord.ipsn.down.labels,
                               nygc_postmortem_spinal_cord.labels, nygc_postmortem_spinal_cord.p53.labels))
nygc_postmortem_spinal_cord.volcano

Functional enrichment terms

nygc_postmortem_spinal_cord.gost <- nygc_postmortem_spinal_cord$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>%  map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.gprofiles <- nygc_postmortem_spinal_cord.gost$result %>% filter(!str_detect(source, "TF|CORUM|HPA|MIRNA")) %>% arrange( -log10(p_value) ) %>% mutate(query = case_when(query == "query_1" ~ " ALS Down ", query == "query_2" ~ " ALS Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.gprofiles_down <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Down ")
nygc_postmortem_spinal_cord.gprofiles_up <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Up ")
down_list <- c(
"synapse",
"axon",
"neuron differentiation",
"cytoskeleton organization",
"neurogenesis",
"transcription by RNA polymerase II",
"regulation of RNA metabolic process",
"protein modification process")
nygc_postmortem_spinal_cord.gprofiles_filt_down <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Down ")  %>% filter(term_name %in% down_list)
up_list <- c(
"regulation of response to DNA damage stimulus",
"Stabilization of p53",
"programmed cell death",
"endoplasmic reticulum",
"lysosome",
"cellular response to stress",
"response to stress",
"translation",
"protein metabolic process",
"mitochondrion")
nygc_postmortem_spinal_cord.gprofiles_filt_up <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Up ")  %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.gprofiles_filt_down, nygc_postmortem_spinal_cord.gprofiles_filt_up) %>% mutate(query = factor(query, levels = c(" ALS Up ", " ALS Down ")))
nygc_postmortem_spinal_cord.gprofiles_filt_combined = nygc_postmortem_spinal_cord.gprofiles_filt_combined %>% mutate(term_name = gsub("regulation of ", "", term_name), term_name = gsub("and","&",term_name), term_name = gsub("cellular ","",term_name), term_name = gsub("response to |transcription by ","",term_name)) 
nygc_postmortem_spinal_cord.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.gprofiles_filt_combined)# + theme_oz() #theme(legend.position = "none", axis.text = element_text(size=10))
nygc_postmortem_spinal_cord.go_curated

Pathways and transcription factors

net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix <- nygc_postmortem_spinal_cord$res %>% select(gene_name, stat) %>%  drop_na(stat) %>% distinct(gene_name, .keep_all = TRUE) %>% tibble_to_matrix(stat, row_names = "gene_name")
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts <- run_wmean(mat=nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = source, group2 = source)
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.pathwayActivity <- ggplot(nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts, aes(x=reorder(source, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + 
    theme_oz() + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    labs(title = "Signaling Pathways", x = "", y = "Normalised enrichment ALS vs CTRL") +
    stat_pvalue_manual(nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts, label = "p.signif", y.position = 6, xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) 
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.pathwayActivity

# p53 weights
nygc_postmortem_spinal_cord_progeny.p53weights = net_progeny %>% filter(source == "p53", target %in% nygc_postmortem_spinal_cord$res$gene_name) %>% left_join(select(nygc_postmortem_spinal_cord$res, target=gene_name,stat)) %>% 
  mutate(color = case_when(weight > 0 & stat > 0 ~ '1', weight > 0 & stat < 0 ~ '2', weight < 0 & stat > 0 ~ '2', weight < 0 & stat < 0 ~ '1', TRUE ~ "3"))
nygc_postmortem_spinal_cord_progeny.p53weights.plot = ggplot(nygc_postmortem_spinal_cord_progeny.p53weights, aes(x = weight, y = stat, color = color)) + geom_point(size = 0.5) +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord_progeny.p53weights, (weight > 10 & stat > 2) | (weight > 0 & stat > 6)), aes(label = target), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.1) +
  theme_oz() + theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
  labs(title = "p53 pathway weights", x = "p53 pathway gene weight", y = "ALS vs CTRL test statistic") +
  geom_vline(xintercept = 0, linetype = 'dotted') +  geom_hline(yintercept = 0, linetype = 'dotted') +
  coord_cartesian(xlim = c(-8,18), ylim = c(-7,10)) 
nygc_postmortem_spinal_cord_progeny.p53weights.plot

# TFs
net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts <- run_wmean(mat=nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 10000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% 
  mutate(sig = case_when(p_value < 0.06 & abs(score) < 4 ~ "sig", p_value < 0.06 & abs(score) >= 4 ~ "sig_strong", p_value > 0.05 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.up = nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts %>% top_n(score, n = 3) %>% pull(source)
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.down = nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts %>% top_n(-score, n = 3) %>% pull(source)

# Volcano of TFs score vs p_value
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.volcano = ggplot(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
  labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "Transcription Factors") +
    guides(colour = "none") +
    theme_oz() +  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
  # scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) + scale_x_continuous(limits = c(-xmax,xmax))
  geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts, source %in% c(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.up, nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.down, "TP53")), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +  
  geom_vline(xintercept = 0, linetype = 'dotted')
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.volcano

Compare postmortem with iPSNs

postmortem_ipsn.als_vs_ctrl = select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_als_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.als_vs_ctrl.overlapping = postmortem_ipsn.als_vs_ctrl %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name) # 50
postmortem_ipsn_de_list = list("Postmortem spinal cord: ALS vs CTRL" = nygc_postmortem_spinal_cord$res, "iPSN: ALS vs CTRL" = ipsc_mn_als_datasets.res)
postmortem_ipsn.als_vs_ctrl.scatter = plot_de_correlation(model_list1 = postmortem_ipsn_de_list, model_list2 = postmortem_ipsn_de_list, mutation1 = "iPSN: ALS vs CTRL", mutation2 = "Postmortem spinal cord: ALS vs CTRL", gene_labels = c(postmortem_ipsn.als_vs_ctrl.overlapping, "CFAP410"), plot_p_value = TRUE)
postmortem_ipsn.als_vs_ctrl.scatter

Compare mutations in postmortem

# heatmap of correlation coefficients
postmortem_mutations_stat = select(nygc_postmortem_spinal_cord.sporadic$res, gene_id, Sporadic = stat) %>% 
  left_join(select(nygc_postmortem_spinal_cord.c9orf72$res, gene_id, C9orf72 = stat)) %>%
  left_join(select(nygc_postmortem_spinal_cord.sod1$res, gene_id, SOD1 = stat)) %>%
  left_join(select(nygc_postmortem_spinal_cord.fus$res, gene_id, FUS = stat)) %>%
  drop_na() %>%
  column_to_rownames("gene_id")
cormat <- round(cor(postmortem_mutations_stat),2)
melted_cormat <- melt(cormat)
get_lower_tri<-function(cormat){ # Get lower triangle of the correlation matrix
  cormat[upper.tri(cormat)] <- NA
  return(cormat)
}
get_upper_tri <- function(cormat){ # Get upper triangle of the correlation matrix
  cormat[lower.tri(cormat)]<- NA
  return(cormat)
  }
upper_tri <- get_upper_tri(cormat)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cormat){ # Use correlation between variables as distance
  dd <- as.dist((1-cormat)/2)
  hc <- hclust(dd)
  cormat <-cormat[hc$order, hc$order]
}
cormat <- reorder_cormat(cormat)
upper_tri <- get_upper_tri(cormat)
melted_cormat <- melt(upper_tri, na.rm = TRUE) # Melt the correlation matrix
postmortem.dge_correlation_ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white") +  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1,1), space = "Lab", name="Pearson\nCorrelation") + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 3)  +
  theme_minimal() + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.ticks = element_blank(),
                          legend.position = "top", legend.direction = "horizontal", legend.title = element_text(size = 8), axis.text=element_text(size=7.5)) + # coord_fixed()+ 
  guides(fill = guide_colorbar(barwidth = 7, barheight = 1, title.position = "top", title.hjust = 0.5))
postmortem.dge_correlation_ggheatmap

p53 pathway and TP53 TF

nygc_postmortem_spinal_cord_mn_mutations_list = list("Sporadic" = select(nygc_postmortem_spinal_cord.sporadic$res,gene_name,stat), "C9orf72" = select(nygc_postmortem_spinal_cord.c9orf72$res,gene_name,stat), "SOD1" = select(nygc_postmortem_spinal_cord.sod1$res,gene_name,stat), "FUS" = select(nygc_postmortem_spinal_cord.fus$res,gene_name,stat))
mutations <- names(nygc_postmortem_spinal_cord_mn_mutations_list)
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{drop_na(.x)}) 
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})

# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord_mn_mutations_list.progeny_scores = map_df(nygc_postmortem_spinal_cord_mn_mutations_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean')
nygc_postmortem_spinal_cord_mn_mutations_list.p53.progeny_scores = nygc_postmortem_spinal_cord_mn_mutations_list.progeny_scores %>% filter(source == "p53") %>% mutate(decoupleR = "p53 pathway") 

# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord_mn_mutations_list.dorothea_scores = map_df(nygc_postmortem_spinal_cord_mn_mutations_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean')
nygc_postmortem_spinal_cord_mn_mutations_list.TP53.dorothea_scores = nygc_postmortem_spinal_cord_mn_mutations_list.dorothea_scores %>% filter(source == "TP53") %>% mutate(decoupleR = "TP53 TF") 

nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53 = bind_rows(nygc_postmortem_spinal_cord_mn_mutations_list.p53.progeny_scores, nygc_postmortem_spinal_cord_mn_mutations_list.TP53.dorothea_scores) %>% 
  mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)

nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53.plot <- ggplot(nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53, aes(x=reorder(mutation, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
  scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + #scale_fill_npg() +
  theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 30, hjust = 1), plot.title = element_text(hjust = 0.5), legend.position = "none") +
  geom_hline(yintercept = 0, linetype = 3) +
  labs(x = "", y = "Normalised enrichment") + #ylim(c(-12,12)) +
  stat_pvalue_manual(nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53, label = "p.signif", y.position = 6, xmin = "mutation", xmax = NULL, size = 4, hide.ns = TRUE) +
  facet_wrap(~decoupleR, scales = "free")
nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53.plot

Table S7 postmortem & ipsn overlap DGE

Overlapping differentially expressed genes in ALS vs CTRL in both postmortem spinal cord and iPSNs

postmortem_ipsc_overlap.res = 
  mutate(select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, padj.postmortem = padj, stat.postmortem = stat, lfc.postmortem = log2FoldChange), direction.postmortem = case_when(padj.postmortem < 0.05 & lfc.postmortem > 0 ~ "up", padj.postmortem < 0.05 & lfc.postmortem < 0 ~ "down", TRUE ~ "non-significant")) %>%
  full_join(mutate(select(ipsc_mn_als_datasets.res, gene_name, gene_id, padj.ipsn = padj, stat.ipsn = stat, lfc.ipsn = log2FoldChange), direction.ipsn = case_when(padj.ipsn < 0.05 & lfc.ipsn > 0 ~ "up", padj.ipsn < 0.05 & lfc.ipsn < 0 ~ "down", TRUE ~ "non-significant"))) %>%
  select(gene_name, gene_id, everything()) %>% 
  filter(padj.ipsn < 0.05, padj.postmortem < 0.05, ((stat.ipsn > 0 & stat.postmortem > 0) | (stat.ipsn < 0 & stat.postmortem < 0)))
kbl(postmortem_ipsc_overlap.res) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), fixed_thead = T) %>% scroll_box(width = "100%", height = "500px")
gene_name gene_id padj.postmortem stat.postmortem lfc.postmortem direction.postmortem padj.ipsn stat.ipsn lfc.ipsn direction.ipsn
ISCU ENSG00000136003 0.0000000 8.311838 0.4136826 up 0.0118787 4.469863 0.1030001 up
RRM2B ENSG00000048392 0.0000000 8.233753 0.7036306 up 0.0167379 4.352172 0.2020891 up
PAQR4 ENSG00000162073 0.0000000 -7.967077 -0.5929364 down 0.0397368 -4.022565 -0.1423603 down
RNASEL ENSG00000135828 0.0000000 6.810334 0.3033616 up 0.0032664 5.089660 0.5962462 up
PTCHD4 ENSG00000244694 0.0000000 6.182531 0.5597501 up 0.0086407 4.663660 0.3292627 up
PLIN5 ENSG00000214456 0.0000001 -5.759645 -0.7603759 down 0.0328980 -4.107068 -0.4650748 down
FBXO22 ENSG00000167196 0.0000067 4.891648 0.2366739 up 0.0461030 3.967703 0.1144354 up
SOD2-OT1 ENSG00000285427 0.0000255 -4.587429 -0.4386346 down 0.0170913 -4.337828 -0.2587170 down
SESN1 ENSG00000080546 0.0001084 4.240309 0.2871041 up 0.0163788 4.377748 0.1711968 up
MTCO2P12 ENSG00000229344 0.0001459 -4.165748 -0.7582162 down 0.0092457 -4.621477 -1.4451766 down
OLIG1 ENSG00000184221 0.0009450 -3.663960 -0.3287542 down 0.0328980 -4.101139 -0.7260811 down
OLIG2 ENSG00000205927 0.0010312 -3.638942 -0.3100081 down 0.0297222 -4.154218 -0.7014318 down
CDKN1A ENSG00000124762 0.0038023 3.243198 0.5623347 up 0.0198326 4.283226 0.3542375 up
CROCC2 ENSG00000226321 0.0088169 -2.961274 -1.0379353 down 0.0253547 -4.206046 -1.5582071 down
RPL29P11 ENSG00000224858 0.0115688 -2.864505 -0.6140480 down 0.0397368 -4.020057 -0.9411368 down
AC025171.1 ENSG00000177738 0.0283717 -2.524708 -0.1924196 down 0.0077084 -4.755944 -0.3937060 down
MTCO1P12 ENSG00000237973 0.0356193 -2.432034 -0.5159713 down 0.0000000 -8.158533 -1.6010349 down

Scatter postmortem & iPSN subgroups

postmortem_ipsc_mn_mutations_list = list("Sporadic iPSN" = ipsc_mn_sporadic_datasets.res, "C9orf72 iPSN" = ipsc_mn_c9orf72_datasets.res, "SOD1 iPSN" = ipsc_mn_sod1_datasets.res, "FUS iPSN" = ipsc_mn_fus_datasets.res, "Sporadic postmortem" = nygc_postmortem_spinal_cord$res, "C9orf72 postmortem" = nygc_postmortem_spinal_cord.c9orf72$res, "SOD1 postmortem" = nygc_postmortem_spinal_cord.sod1$res, "FUS postmortem" = nygc_postmortem_spinal_cord.fus$res)

postmortem_ipsn.sporadic.join = select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_sporadic_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.sporadic.overlapping = postmortem_ipsn.sporadic.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.c9orf72.join = select(nygc_postmortem_spinal_cord.c9orf72$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_c9orf72_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.c9orf72.overlapping = postmortem_ipsn.c9orf72.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.sod1.join = select(nygc_postmortem_spinal_cord.sod1$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_sod1_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.sod1.overlapping = postmortem_ipsn.sod1.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.fus.join = select(nygc_postmortem_spinal_cord.fus$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_fus_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.fus.overlapping = postmortem_ipsn.fus.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)

postmortem_ipsc_mn_mutations_list_t_stat_plot <- plot_grid(
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "Sporadic iPSN", mutation2 = "Sporadic postmortem", gene_labels = postmortem_ipsn.sporadic.overlapping, plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "C9orf72 iPSN", mutation2 = "C9orf72 postmortem", gene_labels = c(postmortem_ipsn.c9orf72.overlapping, "C9orf72"), plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "SOD1 iPSN", mutation2 = "SOD1 postmortem", gene_labels = c(postmortem_ipsn.sod1.overlapping, "SOD1"), plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "FUS iPSN", mutation2 = "FUS postmortem", gene_labels = c(postmortem_ipsn.fus.overlapping, "FUS"), plot_p_value = TRUE),
  ncol = 2, nrow = 2)# & xlim(-8,8) & ylim(-8,8)
postmortem_ipsc_mn_mutations_list_t_stat_plot

Fig S14 Compare mutations postmortem

Compare sporadic, C9orf72, SOD1, FUS

Volcano

postmortem_mutants_bind_datasets.res = bind_rows(mutate(nygc_postmortem_spinal_cord.sporadic$res, mutation = "Sporadic"), mutate(nygc_postmortem_spinal_cord.c9orf72$res, mutation = "C9orf72"), mutate(nygc_postmortem_spinal_cord.sod1$res, mutation = "SOD1"), mutate(nygc_postmortem_spinal_cord.fus$res, mutation = "FUS"))

postmortem_mutants_join_datasets.res = full_join(select(nygc_postmortem_spinal_cord.c9orf72$res, gene_name, padj.c9orf72 = padj, log2FoldChange.c9orf72 = log2FoldChange),
          select(nygc_postmortem_spinal_cord.sporadic$res, gene_name, padj.sals = padj, log2FoldChange.sals = log2FoldChange)) %>%
  full_join(select(nygc_postmortem_spinal_cord.fus$res, gene_name, padj.fus = padj, log2FoldChange.fus = log2FoldChange)) %>%
  full_join(select(nygc_postmortem_spinal_cord.sod1$res, gene_name, padj.sod1 = padj, log2FoldChange.sod1 = log2FoldChange)) #%>% full_join(select(postmortem_vcp_datasets$res, gene_name, padj.vcp = padj, log2FoldChange.vcp = log2FoldChange))
postmortem_mutants_join_datasets.res.labels = postmortem_mutants_join_datasets.res %>% filter(c(padj.c9orf72 < 0.05 & padj.fus < 0.05 & padj.sod1 < 0.05) | 
                                                                                        c(padj.fus < 0.05 & padj.sod1 < 0.05) | 
                                                                                        c(padj.fus < 0.05 & padj.c9orf72 < 0.05)) %>% pull(gene_name)

nygc_postmortem_spinal_cord.fus.labels <- nygc_postmortem_spinal_cord.fus$res %>% filter(-log10(padj)>3, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.fus.volcano <- volcano_plot(nygc_postmortem_spinal_cord.fus$res, "FUS", numerator = "FUS", denomenator = "CTRL", subtitle = "57 Controls vs 2 FUS", ymax = 15, xmax = 4, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.fus.labels))
nygc_postmortem_spinal_cord.fus.volcano

nygc_postmortem_spinal_cord.c9orf72.labels <- nygc_postmortem_spinal_cord.c9orf72$res %>% filter(-log10(padj)>10, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.c9orf72.volcano <- volcano_plot(nygc_postmortem_spinal_cord.c9orf72$res, "C9orf72", numerator = "C9orf72", denomenator = "CTRL", subtitle = "57 Controls vs 36 C9orf72", ymax = 30, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.c9orf72.labels))
nygc_postmortem_spinal_cord.c9orf72.volcano

nygc_postmortem_spinal_cord.sod1.labels <- nygc_postmortem_spinal_cord.sod1$res %>% filter(-log10(padj)>5, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.sod1.volcano <- volcano_plot(nygc_postmortem_spinal_cord.sod1$res, "SOD1", numerator = "SOD1", denomenator = "CTRL", subtitle = "57 Controls vs 5 SOD1", ymax = 17, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.sod1.labels))
nygc_postmortem_spinal_cord.sod1.volcano

nygc_postmortem_spinal_cord.sporadic.labels <- nygc_postmortem_spinal_cord.sporadic$res %>% filter(-log10(padj)>15, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.sporadic.volcano <- volcano_plot(nygc_postmortem_spinal_cord.sporadic$res, "Sporadic", numerator = "Sporadic", denomenator = "CTRL", subtitle = "57 Controls vs 161 Sporadic", ymax = 35, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.sporadic.labels))
nygc_postmortem_spinal_cord.sporadic.volcano

Functional enrichment terms

# c9orf72
nygc_postmortem_spinal_cord.c9orf72.dge_genes <- nygc_postmortem_spinal_cord.c9orf72$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>%  map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.c9orf72.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|MIRNA|HP")) %>% arrange(-log10(p_value)) %>%
  mutate(query = case_when(query == "query_1" ~ " C9orf72 Down ", query == "query_2" ~ " C9orf72 Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_down <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Down ")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_up <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Up ")
down_list <- c(
"regulation of RNA metabolic process",
"ribonucleotide binding",
"cellular localization",
"glutamatergic synapse",
"cytoskeleton organization",
"protein modification process",
"synaptic signaling",
"axon",
"nervous system development",
"synapse",
"protein binding")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"intrinsic apoptotic signaling pathway by p53 class mediator",
"cellular response to DNA damage stimulus",
"Nonsense-Mediated Decay (NMD)",
"apoptotic process",
"programmed cell death",
"cell death",
"mitochondrion",
"cytoplasmic translation",
"response to stress",
"extracellular exosome",
"protein metabolic process",
"protein binding")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Up ")  %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_up) %>%
  mutate(query = factor(query, levels = c(" C9orf72 Up ", " C9orf72 Down ")))
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined$term_name <- gsub("intrinsic apoptotic signaling pathway by p53 class mediator", "apoptotic signaling by p53", nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined$term_name) %>% gsub("regulation of ", "", .) %>% gsub("response to ","",.) %>% gsub(" \\(NMD\\)","",.)
nygc_postmortem_spinal_cord.c9orf72.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.c9orf72.go_curated

# FUS
nygc_postmortem_spinal_cord.fus.dge_genes <- nygc_postmortem_spinal_cord.fus$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>%  map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.fus.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HP|MIRNA")) %>% arrange(-log10(p_value)) %>%
  mutate(query = case_when(query == "query_1" ~ " FUS Down ", query == "query_2" ~ " FUS Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.fus.dge_gprofiles_down <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Down ")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_up <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Up ")
down_list <- c(
"oligodendrocyte differentiation",
"cellular localization",
"protein modification process",
"cytoskeleton organization",
"localization",
"axon",
"cytoskeleton",
"protein binding",
"neuron development",
"synapse",
"nervous system development")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Down ") %>% filter(term_name %in% down_list)
up_list <- c(
# "mitochondrion",
"response to heat",
# "regulation of programmed cell death",
# "response to stress",
"programmed cell death",
"apoptotic process",
"cell death",
"Protein processing in endoplasmic reticulum",
"cellular response to stress",
"response to unfolded protein",
"cytoplasmic translation",
"protein metabolic process",
"protein folding",
# "protein binding",
"extracellular exosome")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Up ")  %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_up) %>%
  mutate(query = factor(query, levels = c(" FUS Up ", " FUS Down ")))
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.)  %>% gsub("regulation of ","",.)
nygc_postmortem_spinal_cord.fus.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.fus.go_curated

# sod1
nygc_postmortem_spinal_cord.sod1.dge_genes <- nygc_postmortem_spinal_cord.sod1$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>%  map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.sod1.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HP|MIRNA")) %>% arrange(-log10(p_value)) %>%
  mutate(query = case_when(query == "query_1" ~ " SOD1 Down ", query == "query_2" ~ " SOD1 Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_down <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Down ")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_up <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Up ")
down_list <- c(
"oligodendrocyte differentiation",
# "cellular localization",
"protein modification process",
"cytoskeleton organization",
# "localization",
"axon",
"cytoskeleton",
"protein binding",
"neuron development",
"synapse",
"nervous system development")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"cellular response to DNA damage stimulus",
# "p53-Independent DNA Damage Response",
"regulation of immune response",
"cytokine production",
"programmed cell death",
"cell death",
"translation",
"endoplasmic reticulum",
"response to stress",
"protein localization",
# "mitochondrial membrane",
"cellular localization",
"lysosome",
"mitochondrion",
# "extracellular exosome",
"protein metabolic process")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Up ")  %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_up) %>%
  mutate(query = factor(query, levels = c(" SOD1 Up ", " SOD1 Down ")))
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.)  %>% gsub("regulation of ","",.) %>% gsub("cellular response to ","",.)
nygc_postmortem_spinal_cord.sod1.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.sod1.go_curated

# sporadic
nygc_postmortem_spinal_cord.sporadic.dge_genes <- nygc_postmortem_spinal_cord.sporadic$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>%  map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.sporadic.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HPA|MIRNA")) %>% arrange(-log10(p_value)) %>%
  mutate(query = case_when(query == "query_1" ~ " Sporadic Down ", query == "query_2" ~ " Sporadic Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_down <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Down ")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_up <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Up ")
down_list <- c(
"histone modification",
"cytoskeleton organization",
"neurogenesis",
"nervous system development",
"organelle organization",
"positive regulation of cellular process",
"protein modification process",
"DNA-templated transcription",
"protein binding")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"cellular response to DNA damage stimulus",
"Stabilization of p53",
# "p53-Dependent G1 DNA Damage Response",
"G1/S DNA Damage Checkpoints",
"Metabolism of RNA",
"regulation of immune system process",
"cell death",
"protein localization",
"endoplasmic reticulum",
"lysosome",
"response to stress",
"cytoplasmic translation",
"protein metabolic process",
"mitochondrion",
"protein binding")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Up ")  %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_up) %>%
  mutate(query = factor(query, levels = c(" Sporadic Up ", " Sporadic Down ")))
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.)  %>% gsub("regulation of ","",.) %>% gsub("cellular response to ","",.)
nygc_postmortem_spinal_cord.sporadic.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.sporadic.go_curated

Upset plot

### Upset plot https://github.com/const-ae/ggupset
postmortem_mutants_bind_datasets.res = bind_rows(mutate(nygc_postmortem_spinal_cord.sporadic$res, mutation = "Sporadic"), mutate(nygc_postmortem_spinal_cord.c9orf72$res, mutation = "C9orf72"),mutate(nygc_postmortem_spinal_cord.sod1$res, mutation = "SOD1"), mutate(nygc_postmortem_spinal_cord.fus$res, mutation = "FUS"))

postmortem_mutants_bind_datasets.res.up.upset = upset_plot(filter(postmortem_mutants_bind_datasets.res, stat > 0), overlap_by = "mutation", overlap_level = "gene_id", split_direction = FALSE, ymax = 2800, order = "freq", order_group_list = c("Sporadic", "C9orf72", "SOD1", "FUS"), title = "Upregulated", label_colours = c("dodgerblue2", "firebrick2", "forestgreen", "gold2"))
postmortem_mutants_bind_datasets.res.up.upset

postmortem_mutants_bind_datasets.res.down.upset = upset_plot(filter(postmortem_mutants_bind_datasets.res, stat < 0), overlap_by = "mutation", overlap_level = "gene_id", split_direction = FALSE, ymax = 2500, order = "freq", order_group_list = c("Sporadic", "C9orf72", "SOD1", "FUS"), title = "Downregulated", label_colours = c("dodgerblue2", "firebrick2", "forestgreen", "gold2"))  
postmortem_mutants_bind_datasets.res.down.upset

Pathways

postmortem_mutations_list = list("Sporadic" = select(nygc_postmortem_spinal_cord.sporadic$res,gene_name,stat), "C9orf72" = select(nygc_postmortem_spinal_cord.c9orf72$res,gene_name,stat), "SOD1" = select(nygc_postmortem_spinal_cord.sod1$res,gene_name,stat), "FUS" = select(nygc_postmortem_spinal_cord.fus$res,gene_name,stat))
mutations <- names(postmortem_mutations_list)
postmortem_mutations_list = map(postmortem_mutations_list, ~{drop_na(.x)}) 
postmortem_mutations_list = map(postmortem_mutations_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
postmortem_mutations_list = map(postmortem_mutations_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})

# net_progeny <- get_progeny(organism = 'human', top = 100)
postmortem_mutations_list.progeny_scores = map_df(postmortem_mutations_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation, source = as.factor(source), source = fct_reorder(source, score), mutation = factor(mutation, levels = c("Sporadic","C9orf72","FUS","SOD1")))

postmortem_mutations_list.progeny_scores.plot <- ggplot(postmortem_mutations_list.progeny_scores, aes(x=mutation, y = score, fill = mutation)) +
  geom_bar(stat='identity', position = "dodge2") +
  scale_fill_npg() +    # scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
  labs(x = "", y = "Signalling Pathways\nALS vs CTRL\nNormalised enrichment") +
  theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90, size = 6), plot.title = element_text(hjust = 0.5), legend.title = element_blank(), legend.position = "none") +
  geom_hline(yintercept = 0, linetype = 3) +
  facet_grid(~source) +
  stat_pvalue_manual(mutate(postmortem_mutations_list.progeny_scores, y.position = case_when(mutation=="Sporadic"~6.5,mutation=="C9orf72"~7,TRUE~7)), label = "p.signif", xmin = "mutation", xmax = NULL, size = 4, hide.ns = TRUE)
postmortem_mutations_list.progeny_scores.plot

Transcription factors

# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
postmortem_mutations_list.dorothea_scores = map_df(postmortem_mutations_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation) %>% 
  mutate(sig = case_when(p_value < 0.1 & abs(score) < 4 ~ "sig", p_value < 0.1 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.1 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction), mutation = factor(mutation, levels = c("Sporadic","C9orf72","SOD1","FUS")))

# Volcano of TFs score vs p_value
postmortem_mutations_merged_dorothea.contrast_acts.volcano = ggplot(postmortem_mutations_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
    labs(y = expression(-log[10]~P~value), x = "Transcription Factors ALS vs CTRL Normalised enrichment") +
    guides(colour = "none") +
    theme_oz() +  #theme(plot.title = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
    geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "C9orf72"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
  geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "FUS"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
  geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "SOD1"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
  geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "Sporadic"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
    geom_text_repel(fontface = "italic", data = filter(postmortem_mutations_list.dorothea_scores,source %in% c("TP53")), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.8) +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  facet_wrap(~mutation, scales = "free", ncol = 4) & ylim(0,2.8)
postmortem_mutations_merged_dorothea.contrast_acts.volcano

Fig S15 TDP43 vs non TDP43 proteinopathy

iPSMN

load(file = here(proj_path, "expression/deseq2/ipsc_mn_tdp43_status.RData"))

ipsc_mn_tdp43_pos_neg_list = list("TDP-43 ALS" = select(ipsc_mn_tdp43_positive_datasets$res,gene_name,stat), "non-TDP-43 ALS" = select(ipsc_mn_tdp43_negative_datasets$res,gene_name,stat))
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{drop_na(.x)}) 
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})

# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
ipsc_mn_tdp43_pos_neg_list.progeny_scores = map_df(ipsc_mn_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)

ipsc_mn_tdp43_positive_datasets_progeny.pathwayActivity <- ggplot(ipsc_mn_tdp43_pos_neg_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + 
    theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    labs(title = "iPSMN signaling", x = "", y = "ALS vs CTRL\nNormalised enrichment") +
    stat_pvalue_manual(mutate(ipsc_mn_tdp43_pos_neg_list.progeny_scores, y.position = case_when(mutation == "non-TDP-43 ALS" & source %in% c("PI3K","TGFb")~ 4, mutation == "TDP-43 ALS" & source %in% c("MAPK") ~ 14, mutation == "TDP-43 ALS"~15, mutation == "non-TDP-43 ALS"~4.5)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
  facet_wrap(~mutation, scales = "free")
ipsc_mn_tdp43_positive_datasets_progeny.pathwayActivity

# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
ipsc_mn_tdp43_pos_neg_list.dorothea_scores = map_df(ipsc_mn_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))

ipsc_mn_tdp43_positive_datasets_dorothea.contrast_acts.volcano = ggplot(ipsc_mn_tdp43_pos_neg_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
  labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "iPSMN TFs") +
    guides(colour = "none") +
    theme_oz() +  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
  geom_text_repel(fontface = "italic", data = filter(ipsc_mn_tdp43_pos_neg_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +  
  geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
  facet_wrap(~mutation)
ipsc_mn_tdp43_positive_datasets_dorothea.contrast_acts.volcano

Postmortem

load(file = here(proj_path, "expression/deseq2/nygc_postmortem_spinal_cord_tdp43_status.RData"))

nygc_postmortem_spinal_cord_tdp43_pos_neg_list = list("TDP-43 ALS" = select(nygc_postmortem_spinal_cord_tdp43_positive_datasets$res,gene_name,stat), "non-TDP-43 ALS" = select(nygc_postmortem_spinal_cord_tdp43_negative_datasets$res,gene_name,stat))
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{drop_na(.x)}) 
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})

# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores = map_df(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)

nygc_postmortem_spinal_cord_tdp43_positive_datasets_progeny.pathwayActivity <- ggplot(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + 
    theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    labs(title = "Postmortem signaling", x = "", y = "ALS vs CTRL\nNormalised enrichment") +
    stat_pvalue_manual(mutate(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores, y.position = case_when(mutation == "non-TDP-43 ALS" & source == "TNFa" ~ 6, mutation == "TDP-43 ALS" & source %in% c("TNFa","NFkB") ~ 5, mutation == "TDP-43 ALS"~5.5, mutation == "non-TDP-43 ALS"~6.5)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
  facet_wrap(~mutation, scales = "free")
nygc_postmortem_spinal_cord_tdp43_positive_datasets_progeny.pathwayActivity

# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores = map_df(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))

nygc_postmortem_spinal_cord_tdp43_positive_datasets_dorothea.contrast_acts.volcano = ggplot(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
  labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "Postmortem TFs") +
    guides(colour = "none") +
    theme_oz() +  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
  geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +  
  geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
  facet_wrap(~mutation, scales = "free")
nygc_postmortem_spinal_cord_tdp43_positive_datasets_dorothea.contrast_acts.volcano

TDP43 knockdown

Liu TDP43+ vs - postmortem TDP43 knockdown models

neuronal_nuclei_tdp43_liu = readRDS(here(proj_path,"expression/deseq2/neuronal_nuclei_tdp43_liu.rds"))
tdp43_kd_datasets = readRDS(here(proj_path,"expression/deseq2/tdp43_kd_datasets.rds"))

ipsc_mn_tdp43_neuron_kd_list = list("Neuronal Nuclei Brain" = select(neuronal_nuclei_tdp43_liu$res,gene_name,stat), 
                                  "TDP-43 knockdown cells" = select(tdp43_kd_datasets$res,gene_name,stat))
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{drop_na(.x)}) 
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})

# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
ipsc_mn_tdp43_neuron_kd_list.progeny_scores = map_df(ipsc_mn_tdp43_neuron_kd_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)

ipsc_mn_tdp43_neuron_kd_datasets_progeny.pathwayActivity <- ggplot(ipsc_mn_tdp43_neuron_kd_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + 
    theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    labs(title = "TDP-43 depletion", x = "", y = "TDP-43 depleted vs TDP-43 retained\nNormalised enrichment") +
     stat_pvalue_manual(mutate(ipsc_mn_tdp43_neuron_kd_list.progeny_scores, y.position = case_when(mutation == "Neuronal Nuclei Brain" & source == "p53" ~ 3.5, mutation == "TDP-43 knockdown cells" & source %in% c("TNFa","VEGF") ~ 15, mutation == "TDP-43 knockdown cells"~16, mutation == "Neuronal Nuclei Brain"~4)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
  facet_wrap(~mutation, scales = "free")
ipsc_mn_tdp43_neuron_kd_datasets_progeny.pathwayActivity

# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
ipsc_mn_tdp43_neuron_kd_list.dorothea_scores = map_df(ipsc_mn_tdp43_neuron_kd_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))

ipsc_mn_tdp43_neuron_kd_datasets_dorothea.contrast_acts.volcano = ggplot(ipsc_mn_tdp43_neuron_kd_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
  labs(y = expression(-log[10]~P~value), x = "Normalised enrichment TDP-43 depleted vs TDP-43 positive", title = "TDP-43 depletion TFs") +
    guides(colour = "none") +
    theme_oz() +  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
  geom_text_repel(fontface = "italic", data = filter(ipsc_mn_tdp43_neuron_kd_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +  
  geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
  facet_wrap(~mutation, scales = "free")
# ipsc_mn_tdp43_neuron_kd_datasets_dorothea.contrast_acts.volcano

TDP43 overexpression

Mouse model from Maor-nof et al

maor_nof_cell_tdp_vs_ctrl = read_csv(here(proj_path, "expression/maor_nof_cell_tdp/GSE162048_D1CON_vs_D1TDP-ExpDiff.csv")) %>% rename(gene_name.mouse = "Gene") %>% left_join(mus.musculus.gene2ens)

maor_nof_cell_tdp_vs_ctrl.genelist = maor_nof_cell_tdp_vs_ctrl %>% select(gene_name.mouse, stat) %>% drop_na() %>% distinct(gene_name.mouse, .keep_all = TRUE) %>% tibble_to_matrix(stat, row_names = "gene_name.mouse")

# PROGENy
net_progeny.mouse <- get_progeny(organism = 'mouse', top = 100)
maor_nof_cell_tdp_vs_ctrl.progeny_scores = run_wmean(maor_nof_cell_tdp_vs_ctrl.genelist, net=net_progeny.mouse, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""))

maor_nof_cell_tdp_vs_ctrl_progeny.pathwayActivity <- ggplot(maor_nof_cell_tdp_vs_ctrl.progeny_scores, aes(x=reorder(source, score), y = score)) +
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + 
    theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    labs(title = "TDP-43 over-expression", x = "", y = "TDP-43 overexpresison vs CTRL\nNormalised enrichment") +
    stat_pvalue_manual(mutate(maor_nof_cell_tdp_vs_ctrl.progeny_scores, y.position = 9.5, group1=source, group2=source), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE)
maor_nof_cell_tdp_vs_ctrl_progeny.pathwayActivity

# Dorothea
net_dorothea.mouse <- get_dorothea(organism='mouse', levels=c('A', 'B', 'C'))
maor_nof_cell_tdp_vs_ctrl.dorothea_scores = run_wmean(maor_nof_cell_tdp_vs_ctrl.genelist, net=net_dorothea.mouse, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))

maor_nof_cell_tdp_vs_ctrl_dorothea.contrast_acts.volcano = ggplot(maor_nof_cell_tdp_vs_ctrl.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
    scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray",  "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
  labs(y = expression(-log[10]~P~value), x = "Normalised enrichment TDP-43 overexpression vs CTRL", title = "TDP-43 overexpression TFs") +
    guides(colour = "none") +
    theme_oz() +  theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
  geom_text_repel(fontface = "italic", data = filter(maor_nof_cell_tdp_vs_ctrl.dorothea_scores, source %in% "Tp53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +  
  geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8))
maor_nof_cell_tdp_vs_ctrl_dorothea.contrast_acts.volcano